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1  Factual  Data 


This  section  contains  a  factual  listing  of  publications,  research  lectures,  amd  students  supported 
in  whole  or  in  part  under  Gintract  No.  AFOS&-89-0167  with  total  funding  support  of  $128,035 
during  the  period  December  15,  1988  to  June  14,  1991.  In  particular,  our  researdi  advances  have 
led  to  15  open  literature  publications,  including  12  in  the  leading  IEEE  and  SIAM  journals.  They 
are  listed  below  and  the  narrative  to  follow  in  Section  2  is  keyed  to  this  list.  Significant  progress 
made  under  the  auspices  of  the  contract  as  well  as  plans  for  the  future  will  be  documented  in  the 
narrative.  Finally,  copies  of  the  abstracts  of  papers  published  or  accepted  for  publication  will  follow 
in  Section  3. 


1.1  Publications  Supported  by  this  AFOSR  Contract 

[1]  Gahinet,  P.,  A.J.  Laub,  C.  Kenney,  and  G.  Hewer,  "Sensitivity  of  the  Stable  Discrete- Time 
Lyapunov  Equation,”  IEEE  Trans.  Aut.  Control,  AC-35(1990),  1209-1217. 

[2]  Pandey,  P.,  C.  Kenney,  and  A.J.  Laub,  "A  Parallel  Algorithm  for  the  Matrix  Sign  Function,” 
Int.  J.  High  Speed  Computing,  2(1990),  181-191. 

[3]  Gahinet,  P.,  and  A.J.  Laub,  “Computable  Bounds  for  the  Sensitivity  of  the  Algebraic  Riccati 
Equation,”  SIAM  J.  Contr.  Opt.,  28(1990),  1461-1480. 

[4]  Pandey,  P.,  C.  Kenney,  A.J.  Laub,  and  A.  Packard,  “Algorithms  for  Computing  the  Optimal 
Hoo  Norm,”  Proc.  29th  IEEE  Conf.  on  Decision  and  Control,  Honolulu,  Hawaii;  December 
1990;  pp.  2628-2629. 

[5]  Pandey,  P.,  C.  Kenney,  and  A.J.  Laub,  “Numerical  Solution  of  Large-Scale  Riccati  Equa¬ 
tions,”  Proc.  Third  Rockwell  Advanced  Control  Systems/Neural  Network/Signal  Processing 
Conf.,  Anahmm,  California;  January  1991;  pp.  100-112. 

[6]  Kenney,  C.,  and  A.J.  Laub,  “Rational  Iterative  Methods  for  the  Matrix  Sign  Function,”  SIAM 
J.  Matrix  Anal.  Appl.,  12(1991),  273-291. 

[7]  Kenney,  C.,  and  A.J.  Laub,  “Polar  Decompontion  and  Matrix  Sign  Function  Condition  Es¬ 
timates,”  SIAM  J.  Sci.  Stat.  Comp.,  12(1991),  488-504. 

[8]  Roy,  S.,  R.H.  Hashemi,  and  A.J.  Laub,  “Square  Root  Parallel  Kalman  Filtering  Using 
Reduced-Order  Local  Filters,”  IEEE  Trans.  Aerosp.  Electr.  Sys.,  27(1991),  276-289. 


[9]  Laub,  A.J.,  “Invariant  Subspace  Methods  for  the  Numerical  Solution  of  Riccati  Equations,” 
in  The  Riccati  Equation,  S.  Bittanti,  A.J.  Laub,  and  J.C.  Willems  (eds.),  Springer-Verlag, 
Berlin,  1991,  pp.  163-196. 

[10]  Pandey,  P.,  C.  Kenney,  A.  Packard,  and  A.J.  Laub,  “A  Gradient  Method  for  Ccunputing  the 
Optimal  j7oo  Norm,”  IEEE  Tmns.  Aut.  Contr.,  AC-36(1991),  887-890. 


[11]  Ghavimi,  A.,  C.  Kenney,  and  A.J.  Laub,  “Local  Convergence  Analysis  ci  Ctmjagtka  Gradient  ,a _ 

Methods  for  Scdving  Algebraic  Riccati  Equations,”  to  appear  in  IEEE  Trmna.  AuL  Contr., 

1992.  _ 

[12]  Gahinet,  P.,  and  A.J.  Laub,  “Algebraic  Riccati  Eqaatfons  and  the  DiMaaee  to  the  Nearest 

Uncontrollable  Pair,”  to  appear  in  SIAM  J.  Contr.  Opt.,  1992.  .  X  |av»jLi~ii57«r^ 
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[13]  Kenney,  C.,  and  A.J.  Laub,  "On  Scaling  Newton’s  Method  for  Polar  Decomposition  and 
the  Matrix  Sign  FVmction,”  to  appear  in  SIAM  J.  Matrix  Anal.  Appl.,  1992  (early  version 
also  appeared  as  "On  Scaling  Newton’s  Method  for  Polar  Decomposition  and  the  Matrix 
Sign  Faction, ”  Proe.  1990  American  Control  Conf.,  San  Diego,  California;  May  1990;  pp. 
2560-2564). 

[14]  Gudmundsson,  T.,  C.  Kenney,  and  A.J.  Lanb,  "Scaling  of  the  Discrete-Time  Algebraic  Riccati 
Equation  to  Enhance  Stability  of  the  Schur  Solution  Method,”  to  appear  in  IEEE  Trans.  Aut. 
Contr.,  1992. 

[15]  Williams,  T.,  and  A.J.  Laub,  "Orthogonal  Canonical  Forms  for  Second-Order  Systems,”  to 
appear  in  IEEE  Trans.  Aut.  Contr. ,  1992  (early  version  also  appeared  as  "Orthogonal 
Canonical  Forms  for  Second-Order  Systems,”  Proc.  1989  American  Control  Conf.,  Pitts¬ 
burgh,  Pennsylvania;  June  1989;  pp.  1621-1622). 

1.2  Major  Invited  Talks  and  Addresses  Supported  by  this  AFOSR  Contract 

1.  State  Space  Computing:  Past,  Present,  and  Future:  Plenary  Lecture  for  the  SIAM  Confer¬ 
ence,  "Control  in  the  Nineties:  Achievements,  Opportunities,  and  Challenges,”  San  FVancisco, 
CA,  May  17-19, 1989;  also  given  as  a  seminar  at  the  Dept,  of  Mechanical  Engineering,  Uni¬ 
versity  of  California,  Irvine,  Apr.  21,  1989. 

2.  The  Matrix  Sign  Function  and  Riccati  Eguationr.  Systems  Research  Center,  University  of 
Maryland,  College  Park,  MD,  Apt.  27, 1989. 

3.  Numerical  Techniques  for  the  Solution  of  Riccati  Equations:  Invited  Plenary  Tutorial  Lec¬ 
ture  for  the  Workshop  on  “The  Riccati  Equation  in  Control,  Systems,  and  Signals,”  Como, 
Italy,  Jun.  26-28, 1989. 

4.  Control  Algorithms  and  Software  Surveg:  Invited  Plenary  Lecture  for  the  3rd  Annual 
Conference  on  Aerospace  Computational  Control,  Oxnard,  California,  Aug.  28-30, 1989. 

5.  Riccati  Equations  and  the  Matrix  Sign  Function:  RIACS,  NASA  Ames,  Moffett  Field,  CA, 
Oct.  17, 1989;  Dept,  of  Electrical  Engineering,  Princeton  University,  Princeton,  NJ,  Oct.  25, 
1989;  Dept,  of  Electrical  Engineering,  University  of  Illinois,  Urbana,  IL,  Nov.  3,  1989. 

6.  Computational  Problems  in  Control  Theory.  Dept,  of  Electrical  Engineering,  University  of 
Pennsylvania,  Philadelphia,  PA,  Oct.  26,  1989. 

7.  The  Matrix  Sign  Function  and  Large-Scale  RUxati  Equations:  Berkeley  Center  for  S)rstems 
and  Contnd,  Spring  Seminar  Series,  University  of  California,  Berkeley,  CA,  May  9,  1990; 
Invited  Plenary  Lecture  for  SIAM  Annual  National  Meeting,  Chicago,  Illinois,  July  18, 
1990. 

8.  IEEE  Control  Systems  Society  Distinguished  Lecture  Series  —  Numerical  Lmear  Algebra 
Problems  tn  Control  Theory.  Ohio  State  Univeiaity,  Columbus,  OH,  October  22, 1990;  Wright 
State  University,  Dayton,  OH,  October  23,  1990. 
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1.3  Graduate  Students  Supported  by  this  AFOSR  Contract 

1.3.1  Ph.D.  Dissertatioiis  Completed 

1.  Pascal  M.  Gahinet 

•  Perturbational  and  Topological  Aspects  of  Sensitivity  in  Control  Theory 

•  December  1989 

•  Presently:  Research  Scientist,  INRIA,  Domaine  de  Voluceau,  Le  Chesnay  (Paris),  France 

2.  Pradeep  Pandey 

•  Numerical  Algorithms  for  Robust  Control  Problems 

•  December  1990 

•  Presently:  Research  Scientist,  Int^ated  Systems,  Inc.,  Santa  Clara,  California 

1.3.2  Other  Graduate  Students  Supported  and  Expected  Completion  Data 

1.  Thorkell  T.  Gudmundsson  (Ph.D.,  Aug.  1992) 

2.  Ali  R.  Ghavimi  (Ph.D.,  Dec.  1992) 

3.  Philip  Papadopoulos  (Ph.D.,  Dec.  1992) 

Citizenship:  Of  the  5  graduate  students  listed  above,  3  are  U.S.  citizens.  Gahinet  (France) 
held  an  F-1  student  visa  as  does  Gudmundsson  (Iceland).  In  addition  to  the  above,  four  other 
Ph.D.  students  also  work  in  the  P.I.’s  research  group:  Thomas  A.  Bryan  (Aug.  1992),  Mark  A. 
Erickson  (Jun.  1993),  John  J.  Hench  (Aug.  1992),  and  Stephen  C.  Stubberud  (Aug.  1992).  All 
are  supported  by  other  contracts  or  fellowships  and  all  are  U.S.  citizens. 

It  should  also  be  mentioned  that  the  P.I.  and  his  students  benefit  enormously  by  being  members 
of  the  Center  for  Control  Engineering  and  Computation  at  UCSB.  This  Center,  under  the  co¬ 
directorship  of  Professor  Petar  Kokotovii  and  the  P.I.,  consists  of  nine  permanent  members  from 
the  Departments  of  Electrical  and  Computer  Engineering,  Chemical  and  Nuclear  Engineering,  and 
Mechanical  and  Environmental  En^neering.  The  cross-departmental  nature  of  the  Center  pves  it 
great  strength  for  its  most  important  task  which  is  to  initiate  and  coordinate  research  projects  rich  in 
opportunities  for  cross-disciplinary  investigations  and  applications  to  industrial,  environmental,  and 
defense  systems.  In  contacts  with  industry,  the  Center  benefits  from  an  unusuaUy  rich  experience 
of  its  members,  covering  a  wide  range  of  technologies. 

2  Narrative 

In  this  section  we  shaO  highlight  research  progress  made  under  this  AFOSR  Contract  and  give  some 
indication  of  current  directions  of  research.  Of  course,  a  more  detailed  description  of  some  of  these 
research  topics  is  contained  in  the  original  proposal. 

The  primary  objective  of  this  project  has  b^n  the  study  of  algorithms  for  large-scale  computa¬ 
tional  problems  aridag  in  control,  filtering,  and  system  theory.  Much  of  our  work  has  concentrated 
on  matrix  Rlccnti  equations  which  are  absdutdy  fundamental  to  the  fidd.  Substantial  progress 
has  been  made  in  other  areas  as  well  and  we  ^ve  the  highlights  s<Hne  trf  the  more  exciting 
contributions  belosr  with  further  narrative  to  follow. 
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(1) :  new  parallel  algorithms  and  successful  implementations  for  several  key  computational  prob¬ 

lems  in  control  and  filtering 

(2) :  significant  breakthroughs  in  the  area  of  error  and  condition  estimation  for  the  algebraic  Riccati 

equation 

(3) :  development  of  an  entire  family  of  iterative  methods,  together  with  a  complete  convergence 

analysis,  for  computing  the  matrix  sign  function;  this  family  of  methods  is  particularly  effec¬ 
tive  for  the  solution  of  large-scale  invariant  subspace  calculations  because  of  its  amenability 
to  implementation  on  parallel  and  vector  computers 

(4) :  the  first  computationally  reliable  method  for  estimating  the  conditioning  of  the  matrix  sign 

function  based  on  the  Frechet  derivative;  this  work  reveals  fascinating  parallels  between  the 
conditioning  of  the  Riccati  equation  and  the  sign  function  of  an  associated  Hamiltonian  matrix 

(5) :  a  new  scaling  strategy  for  Newton’s  method  for  finding  the  sign  of  a  matrix;  this  method 

eliminates  certain  problems  of  determinantal  scaling  while  remaining  nearly  optimal  in  terms 
of  speed 

(6) :  development  and  analysis  of  conjugate  gradient  methods  for  solving  Riccati  equations  and 

general  classes  of  matrix  equations;  necessary  and  sufficient  conditions  for  convergence  have 
been  derived  in  terms  of  the  invertibility  of  the  associated  FVechet  derivative 

(7) :  a  quadratically  convergent  gradient  method  of  determining  optimal  ffoo  norms;  the  method  is 

much  faster  than  current  bisection  methods  and  can  be  extended  to  more  general  perturbation 
problems 

(8) :  enhanced  understanding  of  certain  matrix  "nearness”  problems 

(9) :  enhanced  understanding  of  what  can  and  can  not  be  done  using  reliable  numerical  procedures 

for  matrix  second-order  models 

Each  of  these  results  will  be  described  in  more  detail  below. 

(1) :  We  have  developed  several  new  algorithms  for  parallel  computers  and  enjoyed  successful  im¬ 

plementations  of  some  of  them.  For  example,  we  have  published  research  in  [8]  on  a  number  of 
strategies  and  hierarchies  for  implementing  Kalman  filters  in  a  decentralized  or  parallel  way. 
A  major  contribution  here  is  to  describe  and  antJyze  various  multisensor  network  scenarios 
whose  signal  processing  tasks  are  amenable  to  multiprocessor  implementation.  A  number  of 
extant  strategies  are  unified  and  extended  and  new  algorithms  are  proposed  which  have  the 
potential  for  approximately  linear  speed-up,  are  reasonably  failure-resistant,  and  are  opti¬ 
mized  with  respect  to  communication  bandwidth  and  memory  requirements  at  the  various 
processors  in  the  architecture.  special  feature  of  the  principally  suggested  architecture  is 
the  ability  to  accommodate  parallel  local  filters  of  smaller  state  dimension  than  the  global 
filter.  A  significant  innovation  in  [8]  relative  to  previous  work  is  the  description  of  specific 
implementation  details  for  so-called  square-root  versions  of  filters  in  both  covariance  and  in¬ 
formation  filter  forms.  Another  parallel  algorithm  for  matrix  Riccati  equations  is  discussed 
in  (3)  bdow. 

(2) :  We  have  succeeded  in  deriving  and  extending  computable  error  bounds  for  the  solution  of 

the  algebraic  Riccati  equation  (ARE).  Theee  bounds  are  based  mi  the  themy  devdoped  by 
Kttuiey,  Laub,  and  Wette  {Sya.  Contr.  Lett.,  12(1989),  241-150  for  the  Schnr  method. 
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and  Math,  Contr.  Sig,  Sya.,  3(1990),  211>224  for  Newton  refinement).  These  results  and 
other  algorithms  based  on  invariant  subspaces  are  reviewed  in  an  extensive  survey  paper 
[9]  which  includes  over  230  references.  The  Schur  and  Newton  results  are  complementary 
in  the  sense  that  the  Schur  method  error  bound,  which  is  based  on  an  invariant  subspace 
perturbation  result  of  Stewart  (SIAM  Review^  15(1973),  727-764),  is  needed  to  guarantee  that 
the  computed  solution  is  within  the  region  of  convergence  of  the  Newton  refinement  method.  A 
combination  of  using  the  Schur  method  and  error  bound  as  a  “starter”  for  Newton  refinement 
with  residual  error  estimation,  can  be  easily  implemented  numerically  and  the  necessary 
changes  to  existing  software  are  minimal.  These  results  have  been  extended  to  discrete  and 
singular  Riccati  problems.  For  example,  the  discrete  sensitivity  results  in  [1]  extend  those 
obtained  for  continuous  systems  by  Kenney  and  Hewer  (SLAM  J.  Contr.  Opt.,  28(1990), 
50-69).  Similar  extensions  of  the  Kenney,  Laub,  and  Wette  results  will  appear  in  [14].  Both 
theoretical  and  computable  bounds  are  determined  and  we  note  that  the  discrete-time  case 
turned  out  to  be  somewhat  nontrivial  to  handle.  We  have  also  developed  a  much  deeper 
general  understanding  of  the  sensitivity  of  Riccati  equations.  Based  on  the  Ph.D.  dissertation 
of  Pascal  Gahinet,  we  have  addressed  in  [3]  the  problem  of  determining  computable  boimds 
for  the  condition  or  sensitivity  of  ARE’s.  Specifically,  when  stdving  for  the  unique  symmetric 
nonnegative  definite  solution  of  an  ARE  in  finite  precision  arithmetic,  it  is  crucial  to  know 
topological  properties  of  such  a  solution  when  the  parameter  matrices  of  the  equation  are 
subject  to  perturbation. 

(3) :  Since  the  Schur  method  may  be  impractical  for  the  very  large  Riccati  problems  which  can 

arise,  for  example,  in  distributed  parameter  control  S3rstenu,  we  have  also  made  enormous 
progress  in  extending  the  matrix  sign  function  approach.  This  work  has  its  roots  in  the 
important  extension  of  the  matrix  sign  function  to  generalized  eigenvalue  problems  developed 
by  Gardiner  and  Laub  (Int.  J.  Control,  44(1986),  823-832).  Those  algorithms  are  based  on 
applying  Newton’s  method  to  a  simple  matrix  equation  for  computing  the  sign  of  a  certain 
matrix  and  then  solving  a  certain  linear  system.  This  has  led  further  to  the  search  for 
more  efficient  methods  of  evaluating  the  matrix  sign  function,  and  to  the  development  of  a 
major  paper  [6]  based  on  Fade  approximation  of  a  certain  hypergeometric  function.  This  key 
paper  introduces  a  new  family  of  algorithms,  of  which  the  classical  Newton  iteration  is  but  a 
special  case.  The  algorithms  are  especially  amenable  to  implementation  on  both  parallel  and 
vector  computers  and  a  complete  numerical  analysis,  including  global  convergence  results, 
is  developed  in  [6].  As  part  of  Pradeep  Pandey’s  Ph.D.  dissertation,  vectorized  and  parallel 
versions  of  these  algorithms  have  been  implemented  on  a  Cray  Y-MP  supercomputer  at 
NASA  Ames.  For  reference  we  note  that  even  our  eariy  results  have  shown  that  a  100th- 
order  Riccati  equation  can  be  solved  on  this  machine  in  0.8  sec.  Admittedly,  not  everyone 
has  access  to  a  Cray  —  at  this  moment.  However,  it  is  important  to  bear  in  mind  that  we  will 
soon  have  Cray-type  computing  available  in  desk-top  workstations  in  the  next  few  years  in 
much  the  same  way  we  presently  have  workstations  with  as  much  or  more  computing  power 
than  the  industry-standard  VAX  computer  of  only  a  few  yean  ago.  Descriptions  oi  Cray 
implementations  of  efficient  parallel  partial-fractkm  versions  oi  hif^-mrder  fwmnlas  from  our 
new  family  of  algorithms  have  been  published  in  [2]  and  [5].  In  fact,  in  [5]  we  discuss  the 
numerical  solution  oi  Riccati  equations  of  order  556  (involidng  Hamiltonian  m^ces  of  order 
1112)  in  joint  wwk  with  Rockwell’s  Rocketdyne  Divkikm.  The  problem  derives  from  a  model 
assodated  with  Space  Station  Freedom  in  whidi  278  modes  are  indnded. 

(4) :  Aside  from  matten  ai  just  efficiency,  the  important  aamerical  question  of  the  sensitivity  of 

matrix  sign  solutions  has  also  been  coasidmed.  In  [7]  atefiabie  condition  estimation  procedure 
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is  presented  which  costs  two  extra  sign  function  evaluations.  This  work  is  motivated  by  some 
of  the  fundamental  research  by  Kenney  and  Laub  (SIAM  J.  Matrix  Anal.  AppL,  10(1989), 
191-209)  on  estimating  condition  of  general  matrix-valued  functions.  Future  research  efforts 
will  attempt  to  reduce  the  cost  of  this  condition  estimate,  and  extend  these  finite-dimensional 
Riccati  results  to  infinite-dimensional  operator  problems  arising  from  distributed  parameter 
control  systems. 

(5) :  Also  associated  with  the  matrix  sign  problem  is  the  use  of  scaling  factors  to  accelerate  conver¬ 

gence.  In  analyzing  optimal  scaling  factors  for  the  related  problem  of  accelerating  Newton’s 
method  for  polar  decompositon,  we  have  discovered  that  the  commonly  used  determinantal 
scaling  for  the  sign  function  can  behave  non-optimally  in  some  rather  ordinary  situations.  A 
novel  scaling  strategy  based  on  the  spectral  radius  is  also  flawed  but  the  analysis  in  [13]  shows 
that  the  strengths  of  these  two  procedures  can  be  combined  in  such  a  way  as  to  produce  a 
simple  and  efficiently  realizable  scaling  method  which  is  almost  always  optimal.  This  research 
has,  of  course,  significant  practical  value  for  these  important  matrix  calculations. 

(6) :  We  have  also  made  great  strides  in  developing  and  analyzing  conjugate  gradient  (eg)  methods 

for  solving  Riccati  equations  and  general  classes  of  matrix  equations.  lu  particular,  we  have 
been  able  to  show  that  the  eg  method  converges  in  a  neighborhood  of  a  solution  if  and  only 
if  the  Fr4chet  derivative  of  the  matrix  problem  at  the  solution  is  invertible.  This  meams,  for 
example,  that  the  eg  method  is  convergent  near  the  positive  extremal  solution  of  an  ARE 
because  the  stability  of  the  closed-loop  system  matrix  ensures  that  the  Ft4chet  derivative 
(which  in  this  case  is  just  the  usual  Lyapunov  operator)  is  nonsingular  at  that  solution.  Our  eg 
algorithms  can  be  applied  to  both  symmetric  and  nonsymmetric  Riccati  equations,  including 
those  in  various  “nonstandard”  formats  (e.g.,  certain  additional  terms).  The  methods  can 
also  be  extended  to  a  wide  class  of  general  nonlinear  matrix-valued  equations.  This  is  a  very 
promising  approach  for  very  large-scale  problenu  and  our  first  results  will  be  published  in 
(111- 

(7) :  In  a  related  development,  a  quadratically  convergent  gradient  method  for  finding  optimal  Hoo 

norms  has  been  derived.  Empirical  evidence  shows  that  this  approach  is  much  faster  than 
current  techniques  such  as  bisection  methods,  and  can  easily  be  extended  to  more  general  Hoo 
problems.  These  results  have  been  published  in  [10].  Other  new  technical  characterizations 
of  Riccati  solutions  arising  in  the  problem  appear  in  [4]. 

(8) :  A  key  paper  [12]  will  soon  be  published  in  the  area  of  matrix  “nearness”  problems.  In  this 

work,  a  thorough  mathematical  treatment  is  given  of  the  key  problem  of  determining  the 
nearness  to  uncontroUability  of  a  given  controllable  state-space  model.  The  key  to<d  used 
in  the  analysis  is  a  connection  between  nearness  to  unstabilizability  and  the  behavior  of 
the  unique  symmetric  positive  definite  stabilizing  stdution  of  an  associated  algebraic  Riccati 
equation. 

(9) :  Finally,  an  important  result  relating  to  the  matrix  triples  commonly  found  in  so-called  matrix 

second-order  models  will  appear  in  [15].  The  basic  idea  is  to  establish  which  canonical  forms 
are  obtainable  under  orthogonal  equivalent  for  the  standard  matrix  triffle  consisting  of  a 
mass  matrix,  a  stiflfoess  matrix,  and  a  d*mpi«g  matrix.  Equivalotce  under  orthogonal  trans- 
formathws  is,  of  course,  crucial  for  numerical  rdiabOity.  It  is  established  that  an  arlntrary 
dampiag  model  can  not  be  used  but  that  orthogonal  reduction  of  the  commonly  used  modal 
damidag  modd  can  be  so  reduced. 


7 


Thus  we  coatinue  to  be  excited  about  the  progress  that  has  been  made  and  are  extremely 
enthusiastic  about  the  prospects  and  opportunities  for  further  research. 
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See  attached. 
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Sensitivity  of  the  Stable  Discrete-Time 
Lyapunov  Equation 
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Absiraa—'n*  wMiilvtty  of  the  stable  ditcrttfiiaw  Lyapaao*  cqaa- 
tloa  b  aaalyMtf  Ibroefb  the  tpcciral  aona  of  ibc  iavcna  LyaRoeov 
opemer.  TMa  laada  to  a  tHrecUy  coaiRtHabIc  easyMo-iaicrprci  scasMivtiy 
aieafara.  aa4  aiao  RrevMaa  iasigM  iato  Ibe  coaaeciioR  bctwtea  scasilW- 
tty,  lUbiHy  radiaa,  aad  caabttlaalag  of  iba  cigtoprobleai  of  ibc  opca> 
looR  state  aMlrtk  Tbeae  resalti  an  aa  eateasioa,  to  ibe  Bitcnte-tttae 
can,  of  aaalogoef  nsalta  for  tbe  coatlaaoaa4iaM  LyaRuao*  cquatioa. 

I.  Introduction 


and 


I 


In  (10],  the  spectral  norm  of  the  inverse  Lyapunov  operator, 
given  by 


max 


n^ih 


X0O  II  A";if  +  AXII2 


PROPERTIES  of  the  Lyapuinv  equation  are  frequently  inves¬ 
tigated  through  the  associated  Kronecker  operator.  When 
analyzing  the  sensitivity  of  the  equation,  this  approach  leads  to 
some  measures  of  conditioning  which  are  difficult  to  interpret  in 
terms  of  the  system  parameters.  Moreover,  due  to  the  size  of  the 
Kronecker  product  matrices  (n^  x  where  n  is  the  order  of 
the  matrices  in  the  Lyapunov  equation  itself),  the  evaluation  of 
these  condition  numbm  may  be  problematic  for  large  systems. 

Recendy,  a  new  sensitivity  measure  was  introduced  (i0|  for 
the  suble  continuous-time  Lyapunov  equation 

/4"A'  + (1.1) 

In  (1.1),  the  unknown  and  the  matrices  A  and  H'  are  in 
C’'"  (the  space  of  n  x  n  complex-valued  matrices),  A  is 
suble,  and  denotes  the  Hermitian  transpose  of  A: 
A^  ^  A^.  No  further  assumption  is  made  on  W.  Before  getting 
into  more  detail,  we  need  the  following  definitions.  A  linear 
mapping  0  over  C"”*,  defined  by 

%(X)^A^X+XA  (1.2) 

is  called  a  cominuous-time  Lyapunov  operator.  Throughout  this 
paper,  the  vector  norm  will  be  the  Euclidean  norm,  i.e.,  Ijx^^ 
>>  x^x.  Recall  the  definitions  of  the  Frobenius  norm  and  the 
spectral  norm  (also  called  2-norm)  of  a  matrix  as.  respectively. 

(E  (trace (W«Af))* 

BA/B/-  BAfaB-»«.(Af) 

where  o^{M)  stands  fbr  the  largest  singular  value  of  M. 
Finaf’v,  these  two  nonns  induce  corresponding  norms  on  the 
space  of  linear  operators  over  C"*",  which  are.  respectively. 


is  shown  to  be  a  relevant  measure  of  the  sensitivity  of  (1.1). 
Furthermore,  this  norm  turns  out  to  be  erpial  to  the  spectral 
norm  of  the  solution  H  of  (1.1).  obtained  when  IF  »  /.  It  also 
has  a  simple  interpretation  in  terms  of  the  open-loop  system 
characteristics,  namely,  in  terms  of  the  Lj-norm  of  the  mini¬ 
mally  damped  solution  of  z  *  At.  More  precisely 

II /^Bi  lle'^'^oB^'^f- 

l*ol-i  Jo 

In  this  paper,  we  extend  these  continuous-time  lesultt  to  the 
discrete  case,  namely,  the  discrete-time  Lyapunov^bquation 

Q{X)  =  Q.  where  Q{X)  m  X -xf^XF.  (1.3) 

The  paper  is  organized  as  follows.  In  the  next  section,  we 
review  some  basic  concepts  concerning  (1.3)  and  perform  a 
perturbation  analysis  in  br^  Frobenius  and  spectral  norm.  For 
each  norm,  it  is  shown  that  the  corresponding  norm  of  the 
operator  0~ '  is  a  lelevam  measure  of  the  conditioning  of  (1.3). 
Some  classical  techniques  of  estimating  the  Frobenius  norm  of 
0~'  are  then  piresen^  in  Section  in.  before  introducing  in 
Section  IV  our  new  sensitivity  results,  based  on  the  evaluation  of 
the  spectral  norm  of  Q~'.  Not  only  does  this  norm  have  a 
closed-form  expression  in  terms  of  F.  but  also  it  can  be  calcu¬ 
lated  exactly  by  solving  a  single  Lyapunov  equation.  Section  V 
provides  some  interpretation  of  B  '  I  r  ■"  i^nns  of  the  system 
ntfural  damping,  and  gives  various  bounds  relating  this  norm  to 
characteristics  of  the  matrix  F.  In  ligiN  of  these  results,  we 
discuss  in  Section  VI  which  characteristics  of  F  are  crucial  for 
the  conditioning  of  (1.3).  Finally,  our  conclusions  are  illustrated 
by  a  few  selecfed  numerical  examples  in  Section  VO. 

n.  The  Lyapunov  Equation  por  Discret8-Timb 
Systems 
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For  continuous-lime  systems  deserfeed  by  the  equation  x  - 
Ax  -f  Ar.  the  stability  of  the  matrix  A  can  be  ascertained  via 
(l.l).  Thm  is  an  equivalent  result  for  discreie-tinia  systeriis 
x«4.i  ■  Fxt  +  Cut-  The  correspondiiv  Lyapunov  stability  cri¬ 
terion  involves  what  we  call  the  discrete-time  Lyapunov  opera¬ 
tor.  defined  aa 

a(X)  m  X~  F"XF  (2.1) 

what*  diB  iMfrioaa  F  and  X*  art  in  C"".  In  iha  saquel.  we 
shell  call  any  equniioii  of  the  qrpe 

X-F"XF^Q 
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ABSTRACT 

We  ptopoee  a  new  parallel  algorithm  for  computing  the  sign  function  of  a  matrix. 
The  algorithm  is  based  on  the  Pade  approxunation  of  a  certain  hypergeomet¬ 
ric  function  which  in  turn  leads  to  a  rational  function  approximation  to  the  sign 
function.  Parallelism  is  achieved  by  developing  a  partial  fraction  expansion  nf  the 
rational  function  approximation  since  each  fraction  can  be  evaluated  on  a  sepa¬ 
rate  processor  in  parallel.  For  the  sign  function  the  partial  fraction  expansion  is 
numerically  attractive  since  the  roots  and  the  weights  ace  known  analytically  and 
can  be  computed  very  accurately.  We  also  present  experimental  results  obtained 
on  a  Cray  Y-MP. 

Keywords:  Matrix  sign  function,  parallel  algorithms,  Padd  approximatioa. 


1.  Introduction.  For  a  complex  scalar  z,  with  Jte(z)  ^  0,  the  sign 
function  dpn(z)  is  defined  as: 


sgn{z) 


f  -1-1  when  Re{z)  >  0 
—1  when  iZe(z)  <  0. 


This  definition  can  be  extended  to  matrices  X  whose  eigenvalues  do 

not  lie  on  the  imaginary  axis  in  the  fdlowing  way  [14,4].  Let  X  =  T~^{D  -I- 
N)T,  where  T  is  nonsingular,  D  =diag((ii,...,dp),  and  N  is  nilpotent  and 
commutes  with  D.  Define  the  sign  of  X  by 

8gn{ X)  =  T  diag (  syn(di ),..., sgn{d,)  ]  r~‘ . 


Because  the  columns  of  /  -  sgn{X)  and  l+sgn{X)  form  bases  of  respectivdy 
the  lelt-half-plane  and  right- half-plane  invariant  subspaces  of  X,  the  matrix 


ui 
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Table  2 

Timing  comparison  for  large  test  matrices. 

[10]  P.  Henrk 

Sons,  ! 

[11]  C.  Kenne- 

size  Newton’s  method 

Rational  method 

[12]  A.  J.  Lai 

(sec) 

(sec) 

Auton 

98  1.1423 

0.2675 

[13]  - .  Eff 

298  3.4565 

1.3390 

Contr 
[14]  J.  D.  Rot 

398  8.0332 

3.1150 

by  use 

[15]  G.  Szbgc 

1939. 

can  be  a  problem  if  the  denominator  polynomial  is  ill-conditioned  amd  the 
weights  are  large  and  of  differing  signs.  However,  we  have  adso  shown  that 
for  the  sign  function,  the  roots  and  weight  are  known  analytically  'nd  can 
be  computed  accurately.  Hence,  the  partial  fraction  expansion  is  numeri¬ 
cally  attractive  for  the  sign  function.  Experimental  results  indicate  that  the 
parallel  algorithm  achieves  the  expected  speedup  when  implemented  on  a 
paraUel  machine  like  a  Cray  Y-MP. 
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COMPUTABLE  BOUNDS  FOR  THE  SENSmVITY  OF  THE  ALGEBRAIC 

RICCATI  EQUATION* 

P.  GAHINETt  AND  A.  1.  LAUBt 


AMfacL  Cn  control  or  estimation  theory,  linear-quadratic  optimization  problems  give  rise  to  the 
so-called  matrix  algebraic  Riccati  equation  (ARE).  For  such  problems,  a  crucial  issue  is  the  existence  and 
uniqueness  of  a  symmetric  nonnegative  definite  stabilizing  solution  to  the  ARE,  and  conditions  on  the 
equation  parameters  ate  known  which  guarantee  both.  However,  in  the  context  of  compuutions  in  finite 
precision  arithmetic,  and  with  imperfect  parameter  identification,  it  it  of  concern  whether  the  ARE  retains 
such  a  solution  in  the  proximity  of  a  given  set  of  parameters,  and  how  sensitive  this  solution  is  to  parameter 
variation. 

In  this  paper,  topological  properties,  such  as  openness  of  the  domain  of  existence  and  continuity  with 
respect  to  parameters,  are  established  for  the  symmetric  nonnegative  definite  subilizing  solution.  Compuuble 
sensitivity  estimates  are  also  derived,  which  quantitatively  define  a  region  of  safe  compuution,  in  terms  of 
the  parameters  of  the  equation. 

Key  woeda.  Riccati  equation,  sensitivity,  stabilizability,  compuuble  bounds 

AMSfMOS)  snbiect  classifications.  49E30.  93B3S.  93B40 

1.  Iitrodoctioa.  The  symmetric  algebraic  Riccati  equation  (ARE)  arises  frequently 
in  control  and  estimation  problems.  Consider  the  continuouS'time  ARE  given  by: 

(1.1)  A’'X  +  XA-XFX  +  G  =  0 

where  all  terms  are  matrices  in  R"“"  (real  square  matrices  of  order  n),  and  F  and  G 
are  symmetric,  nonnegative  definite.  The  case  of  complex-vaiued  matrices  is  qualiu- 
lively  similar  to  the  sequel  but  only  the  real-valued  case  will  be  considered  here  since 
it  is  most  commonly  encountered  in  applications.  Under  the  assumption  that  the  pairs 
(A,  F)  and  (G,  A)  are  stabilizable  and  detectable,  respectively,  there  is  a  unique 
nonnegative  definite  symmetric  stabilizing  solution  X  to  (1.1)  (see  [3]  or  [12]).  By  X 
stabilizing  (for  the  pair  (A,  F)),  we  mean  that  A  -  FX  is  stable,  i.e.,  all  its  eigenvalues 
have  striaiy  negative  real  parts. 

Numerical  algorithms  are  now  available  that  solve  the  ARE  efficiently  and  depend¬ 
ably,  provided  the  original  problem  is  sufficiently  well-conditioned  (see  [13]  or  [1]). 
Well-conditioned  means  that  the  solution  X  is  not  greatly  affected  by  small  perturba¬ 
tions  of  the  daU  A,  F,  G.  In  that  case,  and  with  an  appropriate  scaling  of  the  data 
(cf.  [8]),  the  Schur-type  solvere  yield  accurate  solutions  to  (1.1). 

A  natural  question  following  this  preliminary  remark  is  how  to  assess  the  condition¬ 
ing  of  the  symmetric  ARE,  that  is,  its  sensitivity  to  perturbations  of  the  data.  In  other 
words,  if  we  consider  a. perturbed  version  of  (1.1): 

(1.2)  (A-»-AA)^S-t-S(A  +  AA)-S(F+AF)S+G  +  AG-0, 

under  what  conditions  does  (1.2)  keep  a  unique,  nonnegative  definite  stabilizing 
solution  5?  And  can  we  estimate  the  maximum  discrepancy  ||X  -  5||  for  a  given  range 
of  data  perturbations  AA,  AF.  AG? 
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Abstract 

We  piestot  a  gradient  method  for  computing  the  optimal 
norm  Cor  a  general  Ko  control  problem.  This  method  is  much 
fester  than  a  bisection  method  and  the  additional  coat  of  com¬ 
puting  the  gradient  ii  small.  Convargenes  is  predicated  on  the 
emoothnem  of  the  spectral  radine  of  the  product  of  certain  Rk- 
cati  solutiooa.  Hybrid  bieert  ion  •gradient  methods  can  be  used 
in  the  nonamooth  esse. 

1  Introduction 

Considar  the  problem  described  in  [1,  2).  We  are  interested  in 
the  following  linear  system; 

X  3  As  +  SiW  +  Btu, 

M  M  Cis  f>itis  -b  D\2ttt  (1) 

y  *  Cjx  -f  Ojtw  +  Dji«- 


where 

■'»[-;  J]-  ;]■  (« 

The  Hamiltonian  is  deftnad  similarly.  The  following  theorem  is 
essentially  &««  [1], 

Tknomm  1  There  exists  nn  sdminUle  controller  such  that  ||T,«,!|,s  < 
1  iff  the  faUaminf  three  conditions  hold.* 

(i)  Ba  €  domfjtie)  and  Xa  Sie(Sa)  >  0, 
fii^  Ja  €  dam(Rie)  and  y.  Rie(Ja)  >  0, 
ri»;p.:.p(Jf.y.)<o-‘. 

In  the  saqnal  wn  desetibn  thn  behavkr  of  Ba  and  Xa-  By  duality, 
these  comments  also  apply  to  and  1^.  Consular  the  Riceati  equa¬ 
tion  assodatad  with  Ba  in  (2) 

1l(X,  a) :«  AjX.  +  X.A.  -  -I*  C.  -  0.  (4) 


The  problem  is  to  And  a  controUar  which  is  intemally  stabilising  and 
which  minimises  7  ;w  ||T«w||«»,  the  closed-loop  gain  from  to  to  s.  We 
assnme  that  (A,Bi)  and  (A,Bt)  are  stabiUsable,  and  that  (A,C|) 
sad  (A,Ct)  are  deteetabte.  Thms  assumptions  gnarantaa  that  the 
related  Tfs  problem  is  well  deflnad.  Farther,  consider  the  following 
assumptiont: 

Al.  On  =  0,  DUCi  Oisl  «  ( 0  /  ),  D,,  ISf  Of,]  =  (  0  /  ). 

A2.  Dll  -  0- 


These  assumptions  are  from  [1]  where  it  was  noted  that  a  given  system 
could  be  trunsfermed  so  that  it  meats  these  assumptions  via  '‘loop- 
shifting”  [6].  However,  the  tranfermations  to  make  Du  3  0  depend 
on  the  gain  parameter  7.  This  makes  it  harder  to  characterise  the 
behavior  of  the  Riccati  solutions  as  a  feneticn  of  7  and  hence,  is  not 
suitable  fur  the  gradient  mathod  that  we  present.  Therefore,  we  will 
work  directly  with  the  mors  diflcult  ease  in  which  Du  ^  0.  We 
will  refer  to  a  system  that  satislet  aO  the  assumptions  to  be  in  the 
itandsrd  feni-  A  system  fw  which  On  yi  0  will  be  refened  to  aa 
in  the  fcnerel  farm.  For  notational  simpilcity  we  deAne  a  :m  7-*. 
Scalar-  and  matrix-valued  fhnetioas  of  a  will  be  subscripted  with  a. 
If  /  it  such  a  functiaa  then  ff  and  J"  will  denote  its  Arst  and  secMd 
derivative,  raspectivuly,  with  respect  to  a. 

Oh  :■  I  Oil  Ois  1  and  ti  ;■  [  Dj/Ci  1  and  deAne 
the  Hamiltoaiaa  Ba  by 


-i' “(-0, 


‘Suptssted  la  past  by  NsliMsI  Wsace  Iboadoties  (sad  ATOSR)  sades  Oraat 

No.  tem-ttm. 

'Sappastai  la  peM  by  the  Nadsaal  SetaKn  rouedstim  under  Crass  lt» 
PMAISiSSHT  end  the  Ale  fbwn  OUm  ef  Sduellfc  lUmsrch  unda  Ceewmt  Ita, 
AfOdR^MItr. 

'Ssppested  In  pest  by  the  Ceemsl  BecMtc  Cempeey  end  ibe  OC  MIC1U) 
PMgsem  (SMSSIl 


By  aa  aJmieeiUe  solution  wa  will  moan  a  uniqne  poaitivs  sami-deAnite 
stabilising  soiutien.  For  a  a  0  our  aasumptioan  guarantao  e  priori 
that  aa  admissible  Xt  ;•  XiefOe)  2  0  aris^  Bbwuvar,  far  a  >  0  it 
is  not  clear  far  what  vatnaaofa  an  admisaibiasaiatiem  win  exist.  By 
continuity,  Ba  will  be  in  ^am^Rie)  far  some  a  >  0.  Suppose  Ba  and 
Ja  Arst  fsil  to  be  in  demfJtie)  at  a,  and  Op,  respectivriy,  and  deAne 
a.  :3  minfan.a,).  Wa  will  tafsr  to  an  a  €  (0,a.)  as  /eaeMe.  The 
following  theorems,  which  are  proved  in  [5],  chameterim  the  behavior 
of  Ba  and  Xa  for  a  general  pioUam. 

Tbeorwm  2  For  feesMe  a  we  cen  writs  Ga  »  CjCa-  If  V  spans 
the  unehssrpeMs  sahtpoee  e/  the  pair  (A,C|)  then  it  else  spans  the 
unoheervMe  suhspeee  (Aa.Ca).  Ptrther,  k^Xa)  *  V. 

Theorem  t  Let  a.  >  ai  >  as  >  0.  Then  for  e  general  proitem  the 
XieeaU  eohuione  Xa  end  Ya,  and  the  epeetral  radiae  fenetien  pa  are 
centmneus  and  nondeereaeing,  is.,  Xa,  >  Xap,  andpa,  >  pa#- 


2  Gradient  Method 

In  this  seetian  we  reviow  a  gradient  method  (or  computing  7«p(M-  'Be 

will  ssauma  that  condition  (m)  of  Theorem  1  will  bil  before  condition 
(0  or  (it).  This  impliaa  thm  at  tha  optimal  value  p(XaYa)  »  i,  and 
that  Warn  i*  a  root  at  the  equation 

h(a)  ;w  apa  - 1  *  0.  (S) 

Wa  can  And  the  toot  using  tha  Haifoy-sscaal  (or  Newton's)  method. 
The  darivntivua  of  h(a)  an 

V  ■  a/,  and  hf  ■  2/  (6) 

— that  the  dativativas  of  the  spectral  radiM  taction  most. 

Hsiag  thaae  dsialtloae  the  HaBay  sarant  method  yMds 

aiwoiffta,  whan 
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Numerical  Solution  of  Large-Scale  Riccati  Equations 


PRADEEP  PANDEY,  CHARLES  KENNEY,  and  ALAN  J.  LAUB 
Department  of  Etectrieal  and  Computer  Engineering 
Unioeraitg  of  California 
Santa  Barbara,  CA  9S108 

Abstract 

W«  diaeuM  alforithm  for  obtaining  numerical  solutiona  fet  large  ecalc  algebraic  Riccati 
equations.  By  “large-scale”  we  mean  problems  arising  from  models  involving  matrices  of  dimen¬ 
sions  in  the  hundreds  or  perhaps  thousands.  Time  are  numetoos  sources  for  such  problems, 
e.g.,  control  of  large  space  structures,  distributed  parameter  qrstems,  and  interconnected  power 
systems.  Design  of  control  systems  for  such  large  syetems  places  a  signifteani  burden  on  com¬ 
puting  reeoutces  and  may  require  unacceptably  long  computing  time  on  existing  con^uters.  We 
describe  a  parallel  algorithm  for  solving  Riccati  equations  which  is  baaed  on  the  matrix  sign 
function.  We  also  present  numerical  results  obtained  on  a  Cray  y-mp  for  a  structures  model 
currently  under  study  by  Rockwell. 

1  Introduction 

Algebraic  Riccati  equations  (ARE)  play  a  central  role  in  control  theory.  Indeed,  they  are  one 
of  the  most  deeply  studied  nonlinear  matrix  equations  arising  in  mathematics  and  engineering. 
Riccati  equations  atrise  in  a  variety  of  situations  and  th^  role  in  control  and  system  theory  is  wdl 
established. 

Consider  the  fallowing  state-space  model  of  a  linear  system: 

iKO  =  ^  ^ 

where  the  state  x(<)  €  IR**,  the  input  ti(f)  €  and  the  output  y(0  €  IR^.  In  an  LQG  problem 
we  deflne  an  associated  quadratic  cost  function 

y(u) ;» 

Jo 

and  we  want  to  find  a  state-feedback  control  n  such  that  the  closed-loop  system  is  staUe  and  the 
cost  function  J  is  minimized.  This  optimization  probUm  leads  to  the  CoDowing  ARE: 

A^X  +  JCA  -  XBB^X  -h  C^C  m  0,  (2) 

wheroA€  R"’‘",aadC  €  NatnmlasraaptioMorsUbaisnbaityanddelectaba- 

ity  gnninntee  cxistnee  of  a  X^  €  K**’*  which  is  the  naifne  nonnefUiYe  delalte  stnbttiiiag 
solution  of  (2). 

WeaitlntMtedinefldintalforithiBiforsolYiaflnitMeabAREa.  Onrpriasuymolhntionie 

that  in  the  ftttnn  contral  englaeoM  can  inctnesini^  bo  sacpoelod  to  ho  iatsmiod  In  sdviag  pfoUsms 

oflatgidfaneasloao.  At  prmsnt,  problem  sisos  of  twa^  or  osonlBW  "hnndiodi^  can  ho  handled  by 


100 


L. 


Size  of  problem  542 

Time  (seconds)  207 

One-norm  of  residnai  H  3.9093  x  10~* 

Unfortunately,  we  cannot  deduce  accuracy  of  the  computed  solution  from  the  norm  of  the  residual 
alone.  Suppose  5  is  the  true  admissible  solution  for  the  ARE  in  (2).  To  get  a  better  measure  of 
accuracy  we  use  the  following  bound  from  Kenney  et  oL  [14]: 

||J^-S||<2||n-‘||||72||,  (37) 

where  (2  is  the  closed-loop  Lyapunov  operator  defined  by 

Cl(Z)  :=AjZ-i-  ZAe,  where  A-  BB'^X. 

In  [14]  Kenney  et  oL  also  showed  that  ||n~*||  =  ||Z|j  where  Z  €  It'’*''  solves 

AjZ  +  ZAe  +  I  =  0.  (38) 

We  solved  the  above  equation  to  obtain  ||Z||  s  1.362  x  10^.  the  data  we  conclude  that  the 
solution  X  was  computed  to  at  least  6-digit  accuracy. 

6  Summary 

We  have  proposed  a  new  parallel  algorithm  for  solving  lazge-scak  BIccati  equations.  The  algorithm 
is  based  on  a  partial  fraction  expansion  of  a  certain  Pad4  approximation  used  in  computing  the 
matrix  sign  function.  A  partiad  fraction  expansion  allows  a  paralld  implementation  in  which  each 
fraction  can  be  evaluated  on  a  different  processor.  In  geneial,  the  roots  and  weights  of  such  an 
expansion  have  to  be  computed  numerically.  This  can  be  a  problem  if  the  denominator  polyncsnial 
is  ill-conditioned  and  the  weights  are  large  and  of  difihting  signs.  However,  we  have  also  shown  that 
for  the  sign  function,  the  roots  and  weights  are  known  analytically  and  can  be  computed  accuratdy. 
Hence,  the  partial  fraction  expansion  is  numerically  attractive  for  the  sign  function.  Experimental 
results  indicate  that  the  paralld  algorithm  achieves  the  expected  speedup  when  implemented  on  a 
paraUd  machine  like  a  Cray  y-mp. 
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RATIONAL  ITERATIVE  METHODS  FOR 
THE  MATRIX  SIGN  FUNCTION* 

CHARLES  KENNEYt  and  ALAN  J.  LAUBt 


AbdncL  In  this  paper  an  analysis  of  rational  iterations  for  the  matiu  sign  function  is  presented.  This 
analysis  is  based  on  PadI  approximations  of  a  certain  hypergeometrk  fiiitction  and  it  is  shown  that  local 
convergence  results  for  *‘multiplication-nch*'  polynomial  iterations  also  apply  to  these  rational  methods.  Mul- 
tiplication*rich  methods  are  of  particular  interest  for  many  parallei  and  vector  computing  environments.  The 
main  diagonal  Pad!  recursions,  which  include  Newton's  and  Halley's  methods  as  special  ntiwi.  are  globally 
convergent  and  can  be  implemenied  in  a  multiplication-rich  fashion  which  is  computationally  competitive  with 

the  polynomial  recursions  ( which  are  not  globally  convergent).  Other  rational  iteration  schemes  are  also  discuMed. 

including  Laurent  approximations.  Cayley  power  methods,  and  globally  convergent  eigenvalue  assignment 
methods. 


Kay  wards.  Padd  approximation,  matrix  sign  function,  Riccati  equations,  rational  iterations 
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1.  latrodoctkML  It  is  a  classical  result  that  the  algebraic  Riccati  equation  can  be 
solved  by  using  an  invariant  subspace  of  an  associated  Hamiltonian  matrix.  This  motivated 
the  introduction,  by  Roberts  [21]  in  1971,  of  the  matrix  sign  function  as  a  means  of 
finding  the  positive  and  negative  invariant  subspaces  of  any  matrix  X  which  does  not 
have  eigenvalues  on  the  imaginary  axis.  This  and  subsequent  work  [9]  showed  that  the 
matrix  sign  function  could  be  used  to  solve  many  problems  in  control  theory. 

The  sign  of  A*  can  be  defined  constructively  as  the  limit  of  the  Newton  sequence 

(1.1)  =  +  Xo^X, 

(1.2)  sgn(.ir)-  lim  X^. 


Newton’s  method  has  the  pleasant  feature  that  it  is  globally  convergent;  if  X  has  no 
eigenvalues  on  the  imaginary  axis  then  the  limit  in  ( 1 .2 )  exists.  As  a  definition,  however, 
( 1.2)  does  not  reveal  many  of  the  important  properties  of  the  sign  function.  Because  of 
this,  it  is  useful  to  have  an  equivalent  definition  based  on  the  Jordan  canonical  form  of 
X  (see  [4],  [7]).  For  a  complex  scalar  z  with  Re  z  ^  0.  define  the  sign  of  z  by 


(1.3) 


sgn  z 


1  ifRez>0, 
-1  ifRez<0. 


For  a  complex  matrix  X  such  that  .V(  A’)  <=  C”  U  C~  (i.e.,  A*  has  no  eigenvalues  on  the 
imaginary  axis)  let  F  take  X  to  Jordan  form: 


(1.4) 


X’-  r-' 


p 

.0 


0 

s 


r. 


’Raorivtd  By  llM  441100  September  :s.  I9t9;  eceepiid  far  poMigrtoa  (ia  teviirt  liana)  NevBw  15, 
I9W.  TMiiwuitB  eaeigpongd  by  Neuowei  Soeace  FauadMlBa  (ia4  AkFeioiOaeiflfffciiHgc  RimmtB) 
goat  BCaiy«miy,NBiloailScieBMFouiideooojiotDMSi»00H7,ia4  Air  BwotOfceenriiaillcBiioirk 
ooameiAraab«Ml47. 

tDiglBwaiqfgliiiiricrieadCompMierEii»BeefiiigUBiwwlq>efC>Bfcigig.TeBia  CmiBiiii 
9IIOI(lirt4lliirinmiib,iLiti.iJii). 
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POLAR  DECOMPOSITION  AND  MATRIX  SIGN  FUNCTION 
CONDITION  ESTIMATES* 

CHARLES  KENNEYt  and  ALAN  J.  LAUBt 


Abalract.  This  paper  presenis  reliable  condition  estimation  procedures,  based  '  n  Fr^het  derivatives, 
for  polar  decomposition  and  the  matrix  sign  function.  For  polar  decomposition,  the  condition  number  for 
complex  matrices  is  equal  to  the  reciprocal  of  the  smallest  singular  value,  and  rather  surprisingly,  for  real 
matrices  it  is  equal  to  the  reciprocal  of  the  average  of  the  two  smallest  singular  values.  By  using  inverse 
power  methods,  both  of  these  condition  numbers  can  be  evaluated  at  a  fraction  of  the  cost  of  finding  the 
polar  decomposition. 

Except  for  special  cases,  such  as  for  normal  matrices,  the  condition  number  of  the  matrix  sign  function 
does  not  have  such  a  precise  characterization.  However,  accurate  condition  can  be  obtained  by 

using  explicit  forms  of  the  Frtehet  derivative,  or  its  ftnitc.difletetioe  approximation,  with  a  matridal  inverse 
power  method.  These  methods  typically  require  two  extra  sign  function  evaluations,  and  it  is  an  open 
problem  whether  accurate  estimates  can  be  obtained  for  a  fraction  of  a  function  evaluation,  u  is  the  case 
for  the  polar  decomposition.  Related  results  for  the  stable  Lyapunov  equation  and  Nesvton's  method  for 
the  matrix  square  root  problem  are  discussed. 

Key  words,  polar  decomposition,  matrix  sign  function,  conditioning 
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1.  latnrlactioa.  Both  the  polar  decomposition  and  the  matrix  sign  function  play 
important  roles  in  many  matrix  algorithms  [l3-t3],  [8],  [9],  [14],  [17],  The  sensitivity 
of  these  matrix  functions  is  determined  by  the  norms  of  their  Fri^et  derivatives.  More 
specifically,  let  FCA*)  be  a  matrix  function  which  is  continuously  differentiable  at 
X  in  the  sense  that  there  exists  a  linear  matrix  operator  L{’)^L(‘,X,F)  such  that 
for  any  matrix  Z, 


(1.1) 


F{X  +  8Z)~F(X) 
lim - T - 


L{Z). 


Then  L  is  the  FrNhet  derivative  of  F  at  X  and  we  define  the  absolute  and  relative 
condition  numbers  of  F  at  X,  with  respect  to  a  matrix  norm  ||  -  H.  by  X.(F,  X) « ||Li  ■ 
<nax,z„.,llL(Z)||/l|Z||.  and  X,(F.X)-||L111|X||/|1F(X)||  when  IIFlXlIrsQ.  This 
definition  is  consistent  with  the  condition  theory  of  Rice  [16]  (see  [13]).  For  8  small 
and  ||Z||  - 1,  we  see  from  (1.1)  that 

IIAFII  -  ||F(X  +  «Z)-F(X)1|  -  ||L(Z)||«S  «L1|«, 

which  motivates  the  use  of  ||L||  as  a  condition  measure  for  the  mapping  X  -*  F(X). 
As  a  simple  example,  if  F(X)  >  X*  then 


F(X-i-6Z)-F(X) 

8 


ZX+XZ  +  8Z*, 


so  L(Z)-ZX4.XZ  and  ||/.||  S2||X||. 


*  Rcccivud  by  the  cdiion  June  21.  1919;  acccpiud  for  pubWcarioii  (is  rtriisd  feta)  DeesabarA  IMS. 
Thb  raeatch  wm  tupponad  by  National  Sciatica  FouadaUae  (sad  Ak  Fbrea  OBtea  of  SciaaMr  RaMaich) 
graa  ECS  I7.im7,  Nadonal  Saence  Foundatioa  graai  DMSn4«l7,  aad  Air  font  ORea  af  SeiiaMa 
Raaaareh  aoaina  APOSR-I9.0I67. 
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I.  INTRODUCTION 


Square  Root  Parallel  Kalman 
Filtering  Using  Reduced-Order 
Local  RIters 
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Wi  mmmIm  troHoi  Ocrili  111  patolM  KoIomhi  UNcrini 
wllli  iporiol  aM«illMi  itvMi  !•  root 

««tilo«i  !■  kolk  omHmm  MmmUmi  RR«r  hraa. 
ThrMiikwl  Um  ptft  «•  Um  coaMM  wawaioa  of  otiag 

iho  tirw  not”  lo  roll*  to  o  CMnfc;  factoa  A  lyoctal 

fcitamf  thi  MgailoR  arcMIocloio  li  Ife*  OMM;  to  ococaanoM* 


Um  gMil  niia  Tla  atiaMto*  oai  wrirtMra  or  lafotaiaiia 
laairiwo  (**  ikiir  rooti)  ftoai  Umoi  ioRaM4-*r4tr  IWora 
MO  riRatiR  tt  a  oialfal  IRIor  at  *ach  ittR  la  •raoral*  Ua  MMaa 
lloRoRy  oftlaal  wllaaaii  *a<  Hair  aairtam  orror  rorortaaio  or 
Infiraiallia  aBirtao  (or  llatr  naaro  rooti). 
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Deceotnlized  estimation  problemi  have  been  the 
focus  of  great  interest  in  recent  years  in  mulUsensor 
enviionmenis  such  as  integrated  navigation  systems. 
This  stems  from  the  fact  that  using  a  monolithic 
centralized  filter  in  a  multtsensor  system  can  result  in 
severe  computational  problems  due  to  overloading  the 
filter  with  more  data  than  it  can  handle.  Consequently, 
(he  overall  centralize  i  system  may  be  unreliable  and 
suffer  from  poor  accura^  and  subility.  This  leads 
naturally  to  decentralized  processing  configurations. 

Vhrious  decentralized  and  parallel  versions  of  the 
Kalman  filter  have  been  reported  in  the  literature 
{1-10].  Conventional  Kalman  fihering  aigorithan, 
thou^  well  known  to  be  theoretically  ^obally  optimal 
(e.g.,  [11]),  can  also  be  numericaUy  unreliable.  Tb 
remedy  such  problems,  more  numerical  stable  and 
better  conditioned  taplemenutioas  of  the  Kalman 
filtering  algorithma,  such  as  U>D  or  square  root 
formulations,  ate  employed  (12. 13).  Recent  efforts 
have  concentrated  on  decentraUeed  vetsiotts  of  such 
numerically  stable  fOtering  algorithas,  e.|^  (14-18). 

Decentralized  estimation  offers  a<«qeioaa 
advantages  in  many  applications.  Indeed,  in  many 
circumstances  it  provides  the  most  bgic^  feasible 
processiag  schetnesL  For  instance,  in  a  mnltiaensor 
system  in  which  each  individnai  sensor  hu  its  own 
“built-in”  Kabnan  filler,  one  is  interested  in  combining 
the  estimates  from  these  independent  data  sources 
(le.,  the  bnilt-ia  Kalman  filteis)  to  generate  a  global 
estimate  that  will,  ideally,  be  optimaL  nirthermote, 
decentralization  makes  fbr  easy  feult  detectioo  aad 
isolatioa  (3),  since  the  outpnt  of  eadi  local  sensor 
filter  can  be  tested  and,  if  a  sensor  should  fell, 
it  can  be  eiqpeditiouBly  removed  from  the  sensor 
network  before  it  aff^  the  agpegaie  Ster  output 
Also,  decentralizatioo  increases  the  input  data  rates 
significantly  and  yields  moderate  improwasenis  in  the 
throughput 

The  focus  of  several  existing  paralld  Kabnan  fiber 
structures  ((l-^,  6, 7, 9,  K^)  has  been  to  preserve 
the  oveml  gfobel  optimality  of  the  whole  system, 
which  is  definitely  a  desirable  feature  and  senes  as  a 
benchmark  for  other  systems.  While  afi  the  previous 
resttbs  have  amumed  flitt-ortkr  focal  modek  to 
achieve  global  optimaiiiy,  the  intent  of  this  work  is  to 
invettigRta  scenarfos  and  present  atgorkhma  for  wUdt 
the  saaw  can  be  achieved  usiag  rtducai-okkr  focal 
fibers.  In  general,  this  cannot  be  achieved  tor  aibtoary 
gfobal  aad  (lednoed-oidar)  focal  awdelh  mote  it 
said  aboai  this  in  Sacdoa  0.  Abhod0  thb  Mii  das 

fodSd^SaS  nijila  dn  tmH 
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7  Invariant  Subspace  Methods  for  the  Numerical 
Solution  of  Riccad  Equations 

Alan  J.  Laub 


7.1  Introduction 


In  this  tutorial  paper,  an  overview  is  given  ci  psogiess  over  the  past  ten  to  fifteen 
years  towards  reliable  and  effideiu  numerical  soludoo  of  various  typea  of  Riccad 
equations.  Our  anention  will  be  directad  primarily  to  matrix-valued  •tg**— Ric- 

equatioos  and  numerical  methods  for  their  baaed  on  compulinf  baaea 

for  invariant  subspaces  of  certain  associated  Riccad  equadois  arise  in 

modeling  both  continuous-time  and  discrete-dme  sysmais  in  a  wite  vatieqr  of  ap- 
plkadoiu  in  science  and  engineering.  One  can  sonfy  both  algefaraie  eqnmians  and 
difEetendal  or  difference  equadons.  Both  algebnk  and  difibreadal  or  diflStRiion 
equatioos  can  be  further  clusified  accordiag  to  whether  their  coefficient  matrioes 
give  rise  to  so<alkd  symmetric  or  nonsymmetric  equadons.  Symmeaic  Rkcad 
equadons  can  be  further  classified  accordiag  to  whether  or  not  th^  are  definite  or 
indefinite. 

The  rest  of  the  paper  is  organized  as  follows.  A  brief  review  of  “classicar 
methods  is  followed  by  a  summary  ot  the  now-standard  Schur  metfaod,  introduced 
in  1978,  for  solving  algebraic  Riccad  equations.  Eaiensiom  of  the  basic  Schur 
method,  by  means  associated  generalized  eigenvalue  problems,  are  then  de¬ 
scribed  togedier  with  some  applications.  Next,  some  powerful  new  numerical  te- 
sula  relating  to  Riccad  equadons  in  general  wiD  be  described.  These  include  a 
thorough  analysis  of  itetat^  tefinemeu  via  Newton’s  method  (including  a  com¬ 
putable  esdmam  of  the  region  of  convergence),  a  theorem  on  the  teladon  of  enor 
estimates  to  residuals,  enimarion  of  the  condMon  of  algebraic  Riccad  equations, 
and  promising  new  scaling  strategies.  Newton’s  method  fbr  compuring  the  matrix 
sign  function  is  then  described  and  its  implementaiion  fiv  panlM  algorithms  for 
Rkcad  equadons  (on  a  message-passing  hypttcube  compnmr)  is  outlined  This 
method  is  particularly  well  suited  to  paralleiizadan  and  vectotiiadon  and  has  been 
used  suooMsAilly  to  solve  Curly  large  ordar  (several  hnndmd)  problenM.  A  number 
of  generalizations  of  this  bask  iteration  have  exmnded  its  applicabaiiy  to  a  broader 
range  of  problems.  For  example,  generiltadonsof  themsitrixsigafitacdonmdm 
case  of  miutix  pencils  allows  straightfiotwiid  solution  of  (llianm  rime  Rkcad  equar 
dona.  Ftnihetmote,  the  Newton  iteration  itadf  has  been  genwaHind  oonddmtibly 
and  found  to  be  but  a  special  case  of  a  gensial  finally  of  iieradons  for  the  mattfat 
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•  software:  This  is  almost  always,  of  coiine.  the  ultimate  vehicle  of  reliable  tech¬ 
nology  transfier.  Early  attempts  at  a  large  comprohensive  Fbniaii-based  Riccad 
package  (MCMCK;  see  [11])  will  undoubtedly  be  superseded  by  much  moro 
euily  constructed  packages  based  on  software  such  as  matlab  and  its  clones 
and  imiuttors. 

Ackaowtedgenieitta.  Much  of  the  nsemch  surveyed  in  this  paper  represents  col- 
labosadve  eflbra  with  past  and  preaaat  »»«««*■■«««  ifimting  Thrasyvoulos  Pappas, 
Kwee  Hi  Lee.  Bo  Arnold,  Doug  Beodsr.  Judy  Gardiner,  Matt  ^Mae,  Chiu  Choi, 
Pascal  Gahinet,  Pradeep  Pandey,  ThothaU  Oudmundssoii,  and  AU  GhavimL  It  hu 
also  been  a  special  privilege  to  have  die  coUaboraro  with  my  ctd- 

league.  Charles  Ken^,  of  the  .Sciendftc  Cotuputadoo  Laboratory  in  the  Electrical 
and  Computer  Engineeting  Depaitmem  at  UCSB. 

We  are  grateful  for  research  support  fiou  the  Nadooal  Sdeaoe  Roundadon  (and 
AFOSR)  under  Grant  No.  ECS87-18S97.  die  Nadooal  Sdence  Foundation  under 
Grant  So.  DMS88-00817,  mid  the  Air  Foeee  OfBce  of  Scientiilc  Research  under 
Contract  Na  AFOSR-894167. 
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Local  Convergence  Analysis  of  Conjugate  Gradient 
Methods  for  Solving  Algebraic  Riccati  Equations  * 


All  R.  Ghavimi,  Charles  Kenney,  and  Alan  J.  Laub 
Dept,  of  Electrical  and  Computer  Engineering 
University  of  California 
Santa  Barbara,  CA  93106 


Abstract 

Necessary  and  sufficient  conditions  are  given  for  local  convergence  of  the  con¬ 
jugate  gradient  (eg)  method  for  solving  symmetric  and  nonsymmetric  algetouc 
Riccati  equations.  For  these  problems,  the  Ficobenius  norm  of  the  residual  ma¬ 
trix  is  minimised  via  the  eg  method,  and  convergence  in  a  neighborhood  of  the 
solution  is  predicated  on  the  positive  definiteness  of  the  associated  Hessian  ma¬ 
trix.  For  the  nonsynunetric  case,  the  Hessian  eigenvalues  are  determined  by 
the  squares  of  the  singular  values  of  the  closed-loop  Sylvester  operator.  In  the 
symmetric  case,  the  Hessian  eigenvalues  are  closely  related  to  the  squares  of  the 
closed-loop  Lyapunov  singular  values.  In  particular,  the  Hessian  is  positive  defi¬ 
nite  if  and  only  if  the  associated  opmtor  is  nonsingular.  The  invertibility  of  these 
operators  can  be  expressed  as  a  non-cancellation  condition  on  the  eigenvalues  of 
the  closed-loop  matrices.  For  example,  the  stability  of  the  closed-loop  matrix,  for 
the  positive  semi-definite  Riccati  solution,  ensures  the  invertibility  ol  the  Lya¬ 
punov  operator  and  hence  the  convergence  of  the  eg  method  in  a  neighbwhood 
of  that  elution. 


1  Introduction 


When  minimising  a  scalar  function  /  via  the  conjugate  gradient  (eg)  method,  local  conver¬ 
gence  is  equivalent  to  the  Hessian  of  /  being  positive  definite  at  the  point  of  minimisation  [1]. 


’This  fSBsaieh  ess  seppectsd  in  part  by  th*  Natkeal  Sekaes  Ibaadatke  (aad  AFOSR)  uate  Giaat  No 
ECS87-1M97,  by  tbs  Nstfoeal  Seisaes  Fouadatioa  uadsr  Giant  No  DIISMMW617,  aad  by  tbs  Air  Foics 
Ofles  orSeisatMk  IsssaRb  oadst  Coatiaet  No.  AFOSIb«4)167. 
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Algebraic  Riccati  Equations  and  the  Distance  to 
the  Nearest  Uncontrollable  Pair  * 

P.  Gahinet  and  A.J.  Lauh 

Department  of  Electrical  and  Computer  Enpneeiing 
University  of  California,  Santa  Barbara,  CA  93106 


Abstract 

A  connection  is  established  between  nearness  to  unstabilisability  of  a  stabilisabie 
pair  (A,  B)  of  matrices  and  nearness  to  ragularity  of  the  symmetric  positive  definite 
solution  to  an  associated  algebraic  Riccati  equation.  From  this  result,  computable 
upper  and  lower  bounds  are  derived  for  the  distance  of  (A,B)  to  the  nearest  un¬ 
controllable  pair.  Numerical  tests  confirm  the  validity  of  the  method  and  potential 
applications  are  discussed. 


1  Introduction 

When  numerically  assessing  whether  a  pair  of  matrices  {A,B)  €  R"**  x  R"*'  is  con¬ 
trollable  (or  stabilizable),  tests  which  simply  provide  a  yes/no  answer  are  not  entirely 
satisfactory  [17,  18].  Instead,  an  estimate  of  how  far  the  pair  is  from  the  set  of  uncon¬ 
trollable  (respectivdy,  unstabilizable)  pairs  is  more  relevant.  Unfortunatdy,  this  involves 
a  nonconvex  minimization  in  a  space  of  n  dimensions  and  existing  numerical  methods  to 
search  for  minima  often  suffer  from  the  following  limitations: 

•  the  computed  minima  are  only  local, 

•  a  two-dimensional  search  is  necessary  when  complex  perturbations  are  allowed, 

•  the  speed  of  convergence  is  guaranteed  to  be  quadratic  only  in  the  proximity  ci  the 
local  minim*,  and  a  high  computational  overhead  may  thus  be  attached. 

Few  lower  or  upper  bounds  on  the  distance  to  uncontrollabiBty  are  available  in  the 
literature.  Upper  bounds  were  proposed  in  [1]  bat  they  requite  either  forming  the  con¬ 
trollability  matrix,  or  that  A  be  stable.  A  lower  boand  was  obtained  by  Demmd  in  [6]. 

'This  teswich  wu  suppotted  by  the  Nstioaal  Sdice  Poeadatioa  (sad  AFOSR)  uadw  Gnat  No. 
ECSS7-1S89T  sad  the  Air  Force  Office  of  Sdeatific  Resmtch  aadat  Cootiact  No.  AFOSK-SS-OIST.  RepBea 
ehoald  be  sddreesed  to  the  first  author,  curreatly  wiU  WRIA,  Doaedao  do  Volaoasa,  BP  105,  T81S1  Le 
Chansy  Codaa,  Frsaoe. 
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On  Scaling  Newton’s  Method  for  Polar  Decomposition  and  the 

Matrix  Sign  Function  * 

Charles  Kenney  Alan  J.  Laub 

Department  of  Electrical  and  Computer  Engineering 
University  of  California 
Santa  Barbara,  CA  93106 


Abstract 

A  tight  bound  is  given  on  the  speed  of  convergence  of  Newton’s  method  with  optimal  scaling 
for  the  polar  decomposition  of  a  nonsingular  complex  matrix.  Necessary  and  sufficient  conditions 
are  then  derived  which  tell  when  an  approximation  to  the  c^timal  scaling  value  will  give  better 
results  than  the  unsealed  Newton  method.  For  the  related  matrix  sign  problem,  it  is  shown  that 
optimal  scaling  requires  complete  knowledge  of  the  eigenvalues  of  the  original  matrix.  Because 
this  is  impractical,  we  derive  a  family  of  scaling  methods  which  are  optimal  with  respect  to 
partial  eigenvalue  information.  This  family  mcludes  optimal  scaling  as  well  as  a  ‘semHiptimal’ 
scaling  method  based  on  the  dominant  eigenvalues  of  the  matrix  and  its  inverse.  Semi^ptimal 
scaling  can  be  implemented  using  the  power  method  and  gave  nearly  optimal  performance  on 
a  set  of  test  problems.  These  test  problems  also  show  that  a  variety  of  other  commonly  used 
scaling  strategies,  including  spectral  scaling,  determinantal  scaling,  and  2-norm  scaling  can  result 
in  unduly  slow  convergence. 

Ktywordr.  polar  decomposition,  matrix  sign  function,  Newton’s  method,  optimal  scaling. 
AMS(MOS)  subject  classification:  65F35,65P30,15A18. 

Abbreviated  title:  Scaling  Newton’s  method. 


1  Introduction 


The  polar  decomposition  of  a  noosingular  complex  matrix  A  of  order  m  is  a  matrix  pair  (£f,  H)  such  that 
U  is  unitary,  H  is  Hermitian  and  positive  definite,  and  A~UH.  It  A  has  a  singular  value  decomposition, 
A  =  PLQ",  where  P  and  Q  are  unitary  and  E  =  diag{ai,. . .  ,e'm)  with  0  <»«<•••<  Wi, 

U  =  H  =  QEQ".  (1) 


the^ 


However,  it  is  more  efficient  to  compute  the  polar  decompositioD  using  scaled  Newton  recursions  of  the  form. 


^«+l  =  x(7iiAn  +  (7nA^)~^);  Ao^A,  7n  >  0. 
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Scaling  of  the  Discrete-Time  Algebraic  Riccati  Equation  to 
Enhance  Stability  of  the  Schur  Solution  Method  * 

Thorkell  Gudmundaaoa,  Charles  Kenney,  and  Alan  J.  Laub 
Dept,  of  Electrical  and  Computer  Engineering 
University  of  California 
Santa  Barbara,  CA  93106 


Abstract 

A  simple  scaling  procedure  for  discrete>tiine  Riccati  equations  is  introduced.  This  procedure 
instabilities  which  can  be  associated  with  the  linear  equation  solution  step  of  the  gen¬ 
eralised  Schur  method  without  changing  the  condition  of  the  underlying  problem.  A  computable 
bound  for  the  relative  error  of  the  solution  of  the  Riccati  equation  is  also  derived. 


1.  Introduction 

The  Schur  method  [8]  for  solving  discrete-time  algebraic  Riccati  equations  consists  of  transforming  an  as' 
sociated  generalised  eigenvalue  problem  to  real  Schur  form  using  orthogonal  equivalence  transformations, 
followed  by  the  solution  of  a  system  of  linear  equations.  The  orthogonal  transformations  are  numerically 
well-conditioned,  but  recent  work  (12)  has  suggested  that  the  overall  method  can  appear  numerically  unsta¬ 
ble,  even  when  the  original  equation  is  well-conditioned.  This  can  originate  for  two  different  reasons.  One  is 
the  ill-conditioning  of  a  linear  system  of  equations,  and  the  other  is  related  to  scaling  problems  for  the  basis 
vectors  of  a  certain  subspace.  In  this  paper  we  extend  the  work  which  was  done  for  the  continuous-tune 
Riccati  equation  in  (7)  to  the  discrete-time  equation  and  show  that  this  apparent  numerical  instability  can 
be  eliminated  by  a  Kalar  scaling  procedure.  Moreover,  this  analysis  yields  a  good  computable  bound  on  the 
relative  error  of  tbs  solution  of  the  Riccati  problem. 

Our  procedure  is  not  completely  satisfactory,  because  tbs  scalar  involved  is  a  function  of  tbe  solution  to 
be  computed  and  thus  leaves  open  the  questfon  of  how  to  estimate  it  accuratdy  beforehand.  However,  this 
does  not  invalidats  our  main  result  in  any  way,  namely  that  the  Schur  method  is  not  inherently  numerically 
unstable.  In  fact,  the  problem  can  be  circumvented  by  solving  the  unsealed  equation,  using  oor  error  bound 
to  establish  the  aeentaey  of  the  computed  solution,  and  in  ease  that  is  not  satiafiKtocy,  .using  the  inaeeittate 
solution  to  estiinnte  the  optimal  scaling  parameter.  This  estimate  wiD  almost  always  be  siiflldsnUy  close  to 
the  corrset  vahm  to  yMd  an  accurate  result  on  the  second  pass. 
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Abstract 

It  is  shown  in  this  paper  that  a  linear  second-order  system  with  arbitrary  damping  cannot 
be  reduced  to  Hessenberg-triangular  form  by  means  of  orthogonal  transformations.  However, 
it  is  also  shown  that  such  an  orthogonal  reduction  is  always  possible  for  the  modal  damping 
commonly  assumed  for  models  of  flexible  structures.  In  fact,  it  is  shown  that  modally  damped 
models  can  be  orthogonally  reduced  to  a  new  triangular  second-order  Schur  form. 


1  Introduction 


Second-order  models  arise  naturally  in  the  study  of  many  types  of  physical  systems,  with  common 
examples  being  electrical  and  mechanical  networks.  An  application  area  of  great  practical  interest 
for  dynamics  and  control  is  that  of  flexible  space  structures  (FSS)  [2j,  which  are  commonly  repre¬ 
sented  by  second-order  finite  element  models  of  very  high  dimension.  Now,  continuum  models  of 
structures  are,  to  be  sure,  much  more  elegant  (see,  for  example,  [1,  11])  but  it  is  generally  still  the 
case  that  setting  up  the  governing  partial  differential  equations  and  solving  the  resulting  boundary 
value  problems  can  only  be  done  for  relatively  simple  structures.  In  analyzing  a  realistic  structure 
(spacecraft,  airplane,  etc.),  a  continuous  structure  model  is  seldom  feasible  and  common  engineer¬ 
ing  practice  has  been  to  use  some  method  (usually  finite  elements)  to  get  an  approximate  "M  and 
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